* Set up log
cd $ohie
cap log close
global sysdate: disp %tdYYNNDD  date("`c(current_date)'", "DMY")
qui log using 	"./logs/linmte_no_covars_treat_eff_$sysdate.log", replace

* Set up timer
timer clear 1
timer on 1 

/*----------------------------------------------------------------------*/
/* PROGRAM: linmte_no_covars_treat_eff.do				*/
/*									*/
/* PURPOSE:								*/
/* [*]	This code calculates the treated outcomes, untreated outcomes, 	*/
/*	and treatment effects associated with the linear MTE without 	*/
/*	covariates [MTE(p)] and bootstraps these values in order to	*/
/*	obtain the bootstrapped confidence intervals and significance 	*/
/*	stars. The output data is used in the following exhibits of the */
/*	paper (this code does not create the exhibits. It just creates 	*/
/*	the data needed to support the exhibits): 			*/
/* 		- sumstats_Y_num_UOT					*/
/* 		- mte_linear_Y_num					*/
/*									*/
/* OUTPUT:								*/
/* [*]	linmte_no_covars_treat_eff: This output file contains */
/*	all data necessary for replicating the above-mentioned exhibits.*/
/*----------------------------------------------------------------------*/

* Set up display options
clear
set type double
set more off, permanently

*************************************************
* OTHER SETUP			*
*************************************************

* Set up seeds
* This seed stays the same across all codes.
global Y_num_seed = 	6574358

*********************************************
* VARIABLES AND CONTROLS	  	*
*********************************************	

* Input data set	
local indata = "$final/oregonnumhh1.dta"

* Outcomes
local Ys "Y_num"

* Number of bootstrap replications
* The first replication of the bootstrap process has no re-sampling; in other 
* words the first replication always runs on the full original non-bootstrapped 
* sample.
local reps = 1001

* Endogenous variable	
local D		"any_medicaid"

* Instrument
local Z		"Z"

* Output Excel file for bootstrapped statistics
local out_file "linmte_no_covars_treat_eff"

*********************************************
* RUN CODE  	*
*********************************************	

* Begin looping through outcomes
foreach Y of local Ys {
	
	* Initialize an output matrix for storing the results
	matrix treat_eff_`Y' = J(`reps', 100, .)
	
	* Set the seed
	set seed $`Y'_seed
	
	* Initialize a variable to store names for renaming variables
	local renamevars ""

	* Begin bootstrapping loop
	forval rep = 1/`reps' {
		
		* Display to log the outcome and the bootstrap replication
		di "Bootstrap replication #`rep' for `Y'"
		
	*********************************************
	* CALCULATE MTE(p)	  *
	*********************************************
	* We derive the MTE(p) in this code by deriving the ATO(p) and AUO(p) 
	* functions first (similarly to Brinch et al 2015), and then taking 
	* their appropriate derivatives (see paper for details) to obtain the 
	* MTO(p) and MUO(p) functions. We take the difference between MTO(p) 
	* and MUO(p) to get the MTE(p).
	
		* Set up and re-sample
		use "`indata'", clear
			
		* Drop observations where the outcome is missing before 
		* re-sampling for the bootstrapping
		* The reason to drop observations with a missing outcome before 
		* re-sampling is twofold:
		* (a) By dropping missing outcomes, we cannot oversample 
		* individuals with missing values when bootstrapping
		* (b) By dropping missing outcomes, we ensure that every 
		* bootstrap sample is of equal size
		keep if `Y' != . 

		* Re-sample for bootstrapping 
		 quietly if `rep' > 1 {
			bsample
		}
		
		* Initialize iterator for output to matrix 
		local matrixiter = 1
						
		* Compute ATO(p) and AUO(p)
		* We derive ATO(p) and AUO(p) using their closed form 
		* expressions

		* Calculate the BTTO
		* BTTO = E[Y|D=1, Z=0]
		quietly su `Y' if `D' == 1 & `Z' == 0 
		local BTTO = `r(mean)'
		
		* Calculate the ITTO
		* ITTO = E[Y|D=1, Z=1]
		quietly su `Y' if `D' == 1 & `Z' == 1
		local ITTO = `r(mean)'
		
		* Calculate the IUUO
		* IUUO = E[Y|D=0, Z=1]
		quietly su `Y' if `D' == 0 & `Z' == 1 
		local IUUO = `r(mean)'
					
		* Calculate the BUUO
		* BUUO = E[Y|D=0, Z=0]
		quietly su `Y' if `D' == 0 & `Z' == 0 
		local BUUO = `r(mean)'
					
		* Calculate pB
		* pB = P(D=1|Z=0)
		quietly su `D' if `Z'==0
		local pB = `r(mean)'
		
		* Calculate pI
		* pI = P(D=1|Z=1)
		quietly su `D' if `Z'==1
		local pI = `r(mean)'
			
		* Calculate the ATO(p) slope
		local ATO_slope = (`ITTO' - `BTTO')/(`pI' - `pB')
			
		* Calculate the ATO(p) intercept
		local ATO_intercept = `BTTO' - `ATO_slope'*`pB'
			
		* Calculate the AUO(p) slope
		local AUO_slope = (`IUUO' - `BUUO')/(`pI' - `pB')
			
		* Calculate the AUO(p) intercept
		local AUO_intercept = `BUUO' - `AUO_slope'*`pB'
		
		* Calculate Number Treated (N_T)
		quietly count if `D' ==1
		local N_T = `r(N)'
		
		* Calculate Number UnTreated (N_UT)
		quietly count if `D' ==0
		local N_UT = `r(N)'
						
		* Store the ATO(p) and AUO(p) functions in 1x2 matrices
		* The first element of the matrix is the intercept, the second 
		* is the slope	
		
		matrix ATO_matrix = J(1,2,.)
		matrix AUO_matrix = J(1,2,.)
		
		matrix ATO_matrix[1,1] = `ATO_intercept'
		matrix ATO_matrix[1,2] = `ATO_slope'
		matrix AUO_matrix[1,1] = `AUO_intercept'
		matrix AUO_matrix[1,2] = `AUO_slope'
			
		* Compute MTO(p) and MUO(p) using the appropriate .ado files
		
		der_MTO ATO_matrix
		
		der_MUO AUO_matrix
		
		matrix list MTO_matrix
		matrix list MUO_matrix
			
		* Compute the MTE as the difference between MTO(p) and MUO(p)
		matrix MTE_matrix = J(1, 2, .)	
		forval m=1/2 {		
			matrix MTE_matrix[1,`m'] = MTO_matrix[1,`m'] - MUO_matrix[1,`m']				
		}
		
		matrix list MTE_matrix
			
		* Save the MTO(p), MUO(p), MTE(p) slopes and intercepts in the
		* output matrix and assign them names for renaming later
			
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `pB'
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' pB"
		}
		
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `pI'
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' pI"
		}
		
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `=`pI'-`pB''
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' frac_c"
		}
		
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `=1-`pI''
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' frac_nt"
		}
		
		*N_T
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `N_T'
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' N_T"
		}
		
		*N_UT
		matrix treat_eff_`Y'[`rep', `matrixiter'] = `N_UT'
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' N_UT"
		}
		
		* MTO(p) intercept
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MTO_matrix[1,1]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MTO_intercept"
		}
			
		* MTO(p) slope
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MTO_matrix[1,2]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MTO_slope"
		}
			
		* MUO(p) intercept
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MUO_matrix[1,1]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MUO_intercept"
		}
			
		* MUO(p) slope
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MUO_matrix[1,2]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MUO_slope"
		}
				
		* MTE(p) intercept
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MTE_matrix[1,1]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MTE_intercept"
		}
			
		* MTE(p) slope
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MTE_matrix[1,2]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MTE_slope"
		}

		*Extrapolated value
		matrix treat_eff_`Y'[`rep', `matrixiter'] = MTE_matrix[1,1]+0.915*MTE_matrix[1,2]
		local ++matrixiter
		if `rep'==1 {
			local renamevars "`renamevars' MTE_extrapolate"
		}
												
	*********************************************
	* TREATMENT EFFECTS 	*				
	*********************************************	
	* This section uses the treatment_effects.ado file to compute the 
	* treated outcomes, untreated outcomes,and treatment effects		
			
		* Calculate s(pI) = P(Z=1)
		
		quietly su `Z'
		local s_pI = `r(mean)'
			
		* Calculate the treatment effects using treatment_effects.ado 
			
		treatment_effects MTE_matrix MTO_matrix MUO_matrix `pB' `pI' `s_pI'
			
		* Save all of the treated outcomes, untreated outcomes, and 
		* treatment effects for each group of interest (BT, BU, IT, IU, 
		* RIST, RISU, LA, A) in the output matrix and assign them names
		* for renaming later
		
		local groups "BT BU IT IU RIST RISU LA A"
		local outcomes "TO UO TE"
		
		foreach g of local groups {
		
			foreach o of local outcomes {
		
				matrix treat_eff_`Y'[`rep', `matrixiter'] = sc_`g'`o'
				local ++matrixiter
				if `rep'==1 {
					local renamevars "`renamevars' `g'`o'"
				}
			}
		}			
					
		* Decompose the treated outcome of each group into a selection 
		* and treatment effect and save the decompositions in the 
		* output matrix
		* The prefix "sc_" refers to "scalar"
			
		foreach g of local groups {
		
			* Calculation
			
			* Selection effect
			scalar define selection_`g' = sc_`g'UO/sc_`g'TO
			
			* Treatment effect
			scalar define treat_eff_`g' = sc_`g'TE/sc_`g'TO
			
			* Output to matrix and assign names for renaming later
			
			* Selection effect
			matrix treat_eff_`Y'[`rep', `matrixiter'] = selection_`g'
			local ++matrixiter
			if `rep'==1 {
				local renamevars "`renamevars' selection_`g'"
			}
			
			* Treatment effect
			matrix treat_eff_`Y'[`rep', `matrixiter'] = treat_eff_`g'
			local ++matrixiter
			if `rep'==1 {
				local renamevars "`renamevars' treat_eff_`g'"
			}
		}			

		* Calculate the BOLS, IOLS, and RISOLS via differences	

		* BOLS = BTTO - BUUO
		scalar define BOLS = sc_BTTO - sc_BUUO
			
		* IOLS = ITTO - IUUO
		scalar define IOLS = sc_ITTO - sc_IUUO
			
		* RISOLS = RISTTO - RISUUO
		scalar define RISOLS = sc_RISTTO - sc_RISUUO
			
		* Calculate OLS decompositions
		* For each of BOLS, IOLS, and RISOLS we can decompose the OLS 
		* estimate into a selection and treatment effect in two 
		* different ways - using the treated or the untreated within
		* each OLS estimate. This code computes both ways of decomposing 
		* the OLS estimates.
			
		* BOLS decompositions into selection and treatment effects
			
		* BOLS decomposition using BTTE
				
		* Selection effect
		scalar define selection_BOLS_BT = (BOLS-sc_BTTE)/BOLS
		* Treatment effect
		scalar define treat_eff_BOLS_BT = sc_BTTE/BOLS
				
		* BOLS decomposition using BUTE
			
		* Selection effect
		scalar define selection_BOLS_BU = (BOLS-sc_BUTE)/BOLS
		* Treatment effect
		scalar define treat_eff_BOLS_BU = sc_BUTE/BOLS
			
		* IOLS decompositions into selection and treatment effects
			
		* IOLS decomposition using ITTE
			
		* Selection effect
		scalar define selection_IOLS_IT = (IOLS-sc_ITTE)/IOLS
		* Treatment effect
		scalar define treat_eff_IOLS_IT = sc_ITTE/IOLS
		
		* IOLS decomposition using IUTE
			
		* Selection effect
		scalar define selection_IOLS_IU = (IOLS-sc_IUTE)/IOLS
		* Treatment effect
		scalar define treat_eff_IOLS_IU = sc_IUTE/IOLS
			
		* RISOLS decompositions into selection and treatment effects
		
		* RISOLS decomposition using RISTTE
			
		* Selection effect
		scalar define selection_RISOLS_RIST = (RISOLS-sc_RISTTE)/RISOLS
		* Treatment effect
		scalar define treat_eff_RISOLS_RIST = sc_RISTTE/RISOLS
			
		* RISOLS decomposition using RISUTE
		
		* Selection effect
		scalar define selection_RISOLS_RISU = (RISOLS-sc_RISUTE)/RISOLS
		* Treatment effect
		scalar define treat_eff_RISOLS_RISU = sc_RISUTE/RISOLS
		
		* Output BOLS, IOLS, RISOLS and associated decompositions in 
		* the output matrix and assign them names for renaming later
		
			
		local foroutput "BOLS IOLS RISOLS selection_BOLS_BT treat_eff_BOLS_BT selection_BOLS_BU treat_eff_BOLS_BU selection_IOLS_IT treat_eff_IOLS_IT"
		local foroutput	"`foroutput' selection_IOLS_IU treat_eff_IOLS_IU selection_RISOLS_RIST treat_eff_RISOLS_RIST selection_RISOLS_RISU treat_eff_RISOLS_RISU"
		
		foreach f of local foroutput {
		
			matrix treat_eff_`Y'[`rep', `matrixiter'] = `f'
			local ++matrixiter
			if `rep'==1 {
				local renamevars "`renamevars' `f'"
			}
		}
			
	} // end bootstrap loop	
	
	* Save the output matrix
	svmat double treat_eff_`Y'
			
	* Delete extra variables and observations
	keep treat_eff_*
	
	gen obs_num = _n
	drop if obs_num>`reps'
	drop obs_num
				
	* Rename variables
	local i=1		
	foreach var of local renamevars {			
		rename treat_eff_`Y'`i' `var'
		local ++i
	}
	
	* Drop extraneous variables
	drop treat_eff_`Y'*
		
	* Calculate differences of all the treatment effects with the LATO, 
	* LAUO, and LATE
	local teff "TO UO TE"
	
	foreach _te of local teff {
		* Difference
		gen diff_LA`_te'_BT`_te' = 	LA`_te' - BT`_te'
		gen diff_LA`_te'_BU`_te' = 	LA`_te' - BU`_te'
		gen diff_LA`_te'_IT`_te' = 	LA`_te' - IT`_te'
		gen diff_LA`_te'_IU`_te' = 	LA`_te' - IU`_te'
		gen diff_LA`_te'_RIST`_te' = 	LA`_te' - RIST`_te'
		gen diff_LA`_te'_RISU`_te' = 	LA`_te' - RISU`_te'
		gen diff_LA`_te'_A`_te' = 	LA`_te' - A`_te'
		gen diff_BT`_te'_IU`_te' = 	BT`_te' - IU`_te'

	}
		
	* Output the bootstrap statistics
	bootstrap_statistics $output `out_file' `Y'
		
	* Save the bootstrapped file
	save "$final/treat_eff_boot_`Y'.dta", replace
		
} // end outcomes loop
	
timer off 1
timer list 1
local hours = `r(t1)'/3600
di "Computing time is `hours' hours"
		

qui log close
